Onset and universality of drag reduction in the turbulent Kolmogorov flow 
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We investigate the phenomenon of drag reduction in a viscoelastic fluid model of dilute polymer 
solutions. By means of direct numerical simulations of the three-dimensional turbulent Kolmogorov 
flow we show that drag reduction takes place above a critical Reynolds number Rcc- An explicit 
expression for the dependence of Rcc on polymer elasticity and diffusivity is derived. The values 
of the drag coefflcient obtained for different fluid parameters collapse onto a universal curve when 
plotted as a function of the rescaled Reynolds number Re/Rcc- The analysis of the momentum 
budget allows to gain some insight on the physics of drag reduction, and suggests the existence of 
a maximum drag reduction asymptote for this flow. 



When a viscous fluid is kept in motion by some ex- 
ternal driving, a mean flow is established: the ratio be- 
tween the work made by the force and the kinetic en- 
ergy carried by the mean flow is called the drag coeffi- 
cient, or friction factor. This dimensionless number mea- 
sures the power that has to be supplied to the fluid to 
maintain a given throughput. When the flow is lami- 
nar, the drag coefflcient is inversely proportional to the 
Reynolds number. Upon increasing the intensity of the 
applied force the flow eventually becomes turbulent, and 
the drag coefficient becomes approximately independent 
of the Reynolds number [1], therefore substantially larger 
than in the viscous case. 

In 1949 the British chemist Toms reported that the 
turbulent drag could be reduced by up to 80% through 
the addition of minute amounts (few tenths of p. p.m. in 
weight) of long-chain soluble polymers to water. This ob- 
servation triggered an enormous experimental activity to 
characterize this phenomenon (see, e.g., [2, 3, 4, 5, 6]). 
In spite of these efforts, no fully satisfactory theory of 
drag reduction is available yet. However, a recent break- 
through has been the observation of drag reduction in 
numerical simulations of the turbulent channel flow of 
viscoelastic fluids [7] . Most of the features of experimen- 
tal flows of dilute polymer solutions are successfully re- 
produced by these models, even at the quantitative level 
[8]. Despite these advances, the understanding of drag re- 
duction in the experimentally relevant geometry of pipe 
or channel flow is still hindered by the complexity of these 
flows already at the Newtonian level, i.e. in the absence 
of polymers [9] . This consideration motivated us to in- 
vestigate simpler geometries in the hope that this may 
shed some light on the basic physical mechanisms of drag 
reduction (see, e.g., Ref. [10]). 

In this Letter we present the results of an extensive 
numerical investigation of the viscoelastic turbulent Kol- 
mogorov flow. This system has several analogies with the 
turbulent channel flow, while its main distinctive trait 
is the absence of material boundaries. Notwithstanding 



this major difference we will show that drag reduction 
takes place in the Kolmogorov flow as well. Further- 
more, we observe striking quantitative similarities with 
experimental results in wall-bounded flows: this points 
to the conclusion that the basic physical mechanisms of 
drag reduction be substantially independent of the de- 
tailed structure of the flow. 

To describe the dynamics of a dilute polymer solution 
we adopt the linear viscoelastic model (Oldroyd-B) [11] 

dtu+{u-V)u^-Vp + i^oAu+?^^V -(T + F, (1) 

T 

dt(T+{u-'V)(T=(Vu)'^-(7+(T-{Vu)-2^—^+KAcr. (2) 

T 

The velocity field u is incompressible, the symmetric ma- 
trix cr is the conformation tensor of polymer molecules, 
and its trace trcr is a measure of their elongation. The pa- 
rameter T is the (slowest) polymer relaxation time. The 
matrix of velocity gradients is defined as (Vtt)y = diUj 
and 1 is the unit tensor. The solvent viscosity is denoted 
by i^Q and 77 is the zero-shear contribution of polymers to 
the total solution viscosity v = 1^0(1 + rj). The parame- 
ter 77 is proportional to the polymer concentration. The 
diffusive term kAct is added to prevent numerical insta- 
bilities [12]. The constant forcing F = (i^ cos(z/L), 0, 0) 
maintains the system in a statistically stationary state 
characterized by a mean flow (u). Due to the symme- 
tries of F, the only nonzero component of the mean ve- 
locity is (ux)'. it depends on the shear coordinate z alone, 
vanishes at z = ±{t:/2)L, and is even under reflections 
z — z. Its value at z = 0, (ua:)z=Oi will be denoted 
by U. Finally, we establish a short glossary between the 
Kolmogorov flow and the channel flow: F plays the role 
of the pressure gradient, ttL is analogous to the channel 
height, and U is equivalent to the center line velocity. 

In this framework, we have performed a series of nu- 
merical integrations of eqs. (1) and (2) for a set of values 
of forcing intensity F, at flxed both for the Newtonian 
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FIG. 1: Mean velocity profiles for a Newtonian [r] = 0) and a 
viscoelastic simulation (77 = 0.3, El — 0.019) at given forcing 
amplitude F = 1.5. The measured profiles are undistinguish- 
able from (ux) = U cos{z/L) (full lines) in both cases. The 
effect of elasticity is to increase the peak value U with re- 
spect to the Newtonian case: in the present case this cor- 
responds to a reduction of the drag coefficient, defined in 
eq. (3), of about 40%. In the inset, the profiles of the 
Reynolds stress {u^Uz) ~ Ssm{z/L) and the mean polymer 
stress 2t'o7?T~^ ((Ja;z) = —Tsm{z/L). In this case the Reynolds 
stress is reduced upon polymer addition to approximately 70% 
of its Newtonian value, consistently with experimental results 
at comparable drag reduction [14]. The "missing" turbulent 
shear stress is compensated by the contribution of the poly- 
mer stress: the sum of S and T is equal to F in both the 
Newtonian and viscoelastic case. Data result from the nu- 
merical integration of eqs. (1) and (2) in a periodic cube of 
side 27r by means of a fully dealiased pseudospectral code with 
64^ collocation points. The mean flow lengthscale is L = 1 
and the viscosity is u — 0.015625. Starting from an initial 
configuration with a small amount of energy on the smallest 
modes, after the system evolved into a statistically station- 
ary state, time averages over 100 to 1000 eddy-turnover times 
have been performed to obtain the mean velocity profiles. 



and the viscoelastic case. Comparing results at a given 
F is equivalent to keeping an imposed pressure gradi- 
ent - therefore a fixed wall-shear stress - in channel flow 
experiments (see, e.g. Ref. [13]). We have measured 
the mean profiles of several relevant observables, includ- 
ing the average velocity (u^), the turbulent shear stress 
(Reynolds stress) {u^Uz), and the mean polymer stress 
2h'QriT~^{axz)- The mean flow is accurately described 
by the sinusoidal profile {ux) = U cos(z/L), both in the 
Newtonian and in the viscoelastic fiow [15]. However, as 
shown in Fig. 1, in the latter case the centerline veloc- 
ity U is definitely larger: this is the hallmark of drag 
reduction. It has to be remarked that - at variance 
with wall-bounded flows where drag reduction is always 
accompanied by a structural change in the proflle (see 
e.g. Ref. [3]) - in the Kolmogorov flow the increase in 
throughput takes place just by means of an overall rescal- 
ing of the mean velocity. This is due to the different 



boundary conditions: in channel flows, the proflle in the 
viscous sublayer is left unchanged upon polymer addition 
while the bulk flow increases substantially. This requires 
a reshaping of the mean proflle, that takes actually place 
through the increase of the extent of the buffer region 
(see e.g. Ref. [13]). In the Kolmogorov flow there is no 
constraint on velocity proflles, and drag reduction does 
not necessarily entail their structural change. 

To quantify the effect of viscoelasticity on the mean 
flow, we have deflned the drag coefflcient as 



/ = - 



(3) 



and measured its dependence on the Reynolds number 
Re = UL/v [16]. For Re < \f2 the flow is laminar 
with mean velocity U — FL"^ /v, giving a drag coeffl- 
cient f = Re~^. At Re > 50 the system is already in a 
fully developed turbulent state. For a Newtonian fluid, 
numerical data show that the drag coefflcient is approx- 
imately independent of Re (see Fig. 2). This behaviour 
agrees with the following classical Kolmogorov argument: 
since the average energy input e = FU/2 scales as 
£ = § /L in fully developed turbulence, eq. (3) yields 
a constant drag coefflcient f — (3. The Newtonian mo- 
mentum budget gives F^ = dz{uxUz) (the viscous contri- 
bution being negligible) and therefore a Reynolds stress 
{uxUz) = Ssm{z/L) with S = /3C/^. For the turbulent 
Kolmogorov flow, /3 ~ 0.15. 

When polymers are added / may be reduced with re- 
spect to its Newtonian value, depending on the polymer 
elasticity El = vrjl?^ the Schmidt number Sc — v/k, 
and the concentration r/, as shown in Fig. 2. For the 
highest Reynolds number we can attain in our simula- 
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FIG. 2: The drag coefficient for different viscoelastic fluid 
parameters. Data have been collected from numerical simu- 
lations at different forcing amplitudes F and viscoelastic pa- 
rameters 77, r, K. The statistical error in the determination of 
/ and Re is of the order of the symbols' size. 
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tions the friction factor is reduced by 75%. Drag reduc- 
tion is accompanied by changes in the velocity field sim- 
ilar to those occurring in channel flow experiments and 
simulations: the level of transverse fluctuations {uD is 
reduced while longitudinal fluctuations {(ux — {ux)Y) in- 
crease and high streamwise velocity streaks are observed 
(see Fig. 3). Incidentally, we notice that drag reduction 
is observed at Reynolds numbers definitely smaller than 
the typical experimental values: this is possible thanks 
to the relatively high value of elasticity utilized in our 
simulations. Comparable parameters have been used in 
numerical simulations of the channel flow as well (see, e.g. 
Ref. [7]), and produced a similar effect on the threshold 
for drag reduction. 

From the inspection of Fig. 2 we notice some system- 
atic trend: at moderate Reynolds numbers {Re < 200) 
viscoelastic effects do not alter substantially the value of 
the drag coefficient; at larger Re polymers with a higher 
elasticity are more effective as drag-reducing agents; con- 
versely, polymers with higher diffusivity are less effec- 
tive. To understand the variation of the drag coefficient 
with fluid parameters, we sought a dependence of the 
form / = ip(Re/ Rec) where Rec{El, Sc^r]) is the crit- 
ical Reynolds number for the onset of drag reduction. 
To obtain an explicit expression for Rec we need to ex- 
tend the argument given in Ref. [17] to the case of finite 
polymer diffusivity. The reasoning goes as follows: for 
polymers to be substantially elongated, stretching must 
prevail over elastic relaxation and diffusivity [18]; at the 
onset, the terms appearing in cq. (2) must then satisfy 
(Vu)c ~ 2/t -|- k/L'^] since the transition is incipient we 
can estimate the typical velocity gradient as {ec/i^Y^^, 
and utilizing the expression ec oa / L we finally obtain 

For vanishing diffusivity we recover the result of Ref. [17]. 
Extracting the explicit dependence on polymer concen- 




FIG. 3: Snapshots of the isosurfaces Ux{x,y,z) — U for a 
Newtonian (left) and a viscoelastic simulation (right. El = 
0.019, Sc = 0.016, rj = 0.5). The Reynolds number is 
Re « 350. The arrows show the direction of the mean flow. 
Small-scale turbulent fluctuations, responsible for kinetic en- 
ergy dissipation, are suppressed in the viscoelastic case. A 
high-speed streak in the form a streamwise oriented tube is 
visible in the viscoelastic case (right). 



tration, we have Rcc oc (1 -I- r])~'^/^ , which is compatible 
with the weak dependence on concentration found in ex- 
periments [3]. 

In Fig. 4 we present the same data as in Fig. 2, now 
plotted against the rescaled Reynolds number Re/ Rec- 
The good quality of the collapse supports the validity of 
the relation / = ip{Re/ Rec). The function ip is universal 
with respect to the choice of fluid parameters. Its shape 
will be derived in the following, with the aid of simple 
assumptions, starting from the equation for momentum 
conservation (see Ref. [19] for a similar approach to wall- 
bounded flows). 

Upon time averaging, eq. (1) reduces to — 
~vodl{ux) + dz{{uxU^) - 2iyor]{axz)/T). Utilizing the 
numerical observation that the Reynolds stress S = 
(uxUz) ~ Ssin{z/L) and the polymer stress T = 
2i>Qr]{axz) /t = —Tsm{z/L), we obtain the momentum 
budget F = vaU/L^ + S/L + T/L. The contribution 
lyoU/L^ is relevant only in the laminar regime, and can 
therefore be neglected. The dependence of the stresses 
on the rescaled Reynolds number is presented in Fig. 5. 
Below the threshold the polymer stress is vanishingly 
small whereas the Reynolds stress is 5* ~ /3J7^ ~ FL 
in agreement with the observation of a i?e-independent 
drag coefficient. Above Rcc, the polymer stress makes 
a significant contribution to the momentum budget. At 
the largest Re we can attain, the elastic stress reaches 
almost 50% of the total stress, not far from experimen- 
tal results [20]. Rescaling the stresses with the critical 
velocity squared shows that above the onset 5* tends to 
a constant value 7//^ (see Fig. 5(b)), and the polymer 
stress follows the law T = 5U^ (Fig. 5(c)). The physi- 
cal interpretation of these observations is that above the 
onset of drag reduction an increasing fraction of the mo- 
mentum injected by the external force is sequestered by 
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FIG. 4: The drag coefRcient plotted as a function of the 
rescaled Reynolds number Re/Rec- Symbols as in fig. 2. The 
full line is eq. (5) with /3 = 0.15, 7 = 0.2 and 5 = 0.02. 
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greater ease to velocity fluctuations than to elastic ones. 
Understanding the reasons for this difference would dis- 
close the basic physical mechanisms of drag reduction. 

Remarkably, eq. (5) predicts a maximum drag reduc- 
tion asymptote (see, e.g., Refs. [3, 6]). Indeed, by increas- 
ing the concentration, drag cannot be reduced below the 
asymptote f = S, independently of polymer elasticity and 
diffusivity. In this ultimate regime momentum transfer 
would take place only through polymer stresses. How- 
ever, the present data do not cover a sufficient span of 
values of Re to allow us to confirm or reject this predic- 
tion. Numerical simulations at higher resolution should 
allow to settle this issue. 
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the support of EU under the contract HPRN-CT-2002- 
00300 and of miur-cofin 2001023848. Numerical simula- 
tions have been performed at CINECA (INFM parallel 
computing initiative). 



FIG. 5: (a): Peak values of the Reynolds stress {u^Uz) = 
Ssm{z/L), and of the polymer stress 2r]UoT~^ (a^z) = 
— rsin(2;/L), nondimensionalized by the total stress FL. The 
sum {S + T) /(FL) nearly equals unity for each couple of dat- 
apoints, confirming that the viscous (solvent) stress vqU/L^ 
is negligible at the present Reynolds numbers. (&) The 
Reynolds stress S and (c) the polymer stress T nondimen- 
sionalized by the squared critical velocity Uc'- The full lines 
are: S/U'^ = 7 (left, horizontal), S/U'^ = PiRe/Re^f (left, 
oblique), T/U"^ = S{Re/Reof (right). The numerical param- 
eters are /3 = 0.15, 7 = 0.2 and 5 = 0.02. 



polymers, which are however less effective in absorbing 
it than transverse velocity fluctuations {6 < f3). This re- 
sults in an enhancement of the mean flow with respect 
to the Newtonian case, i.e. drag-reduction. Inserting the 
empirical expressions for S and T, the momentum bud- 
get above the onset reads F = ^U'^/L + SU"^ / and the 
resulting drag coefhcient is 



/ = 



7 



Rec 
Re 



for Re < Rec , 
6 for Re > Rec • 



(5) 



This expression is compared with numerical results in 
Fig. 4, where the values of the parameters 7 and S have 
been obtained from the data shown in Fig. 5. The agree- 
ment is excellent, except possibly for Re « Rec, where 
eq. (5) predicts an abrupt transition: from Fig. 5 this 
rather appears to be a smooth crossover, whose actual 
shape cannot be extracted by means of simple arguments. 
The actual values of /3, 7 and S are not of utmost impor- 
tance since they are likely to depend on the details of the 
driving force, and therefore on the shape of the velocity 
profile. What is crucial to drag reduction is that 6 < (3, 
or - in plain words - that momentum is transferred with 
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